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Abstract 

Gibbs  oscillation  can  show  up  near  flow  regions  with  strong  temperature  gradients  in  the 
numerical  simulation  of  nonhydrostatic  (NH)  mesoscale  atmospheric  flows  when  using  the  high- 
order  discontinuous  Galerkin  (DG)  method.  We  propose  to  incorporate  localized  Laplacian 
artificial  viscosity  in  the  DG  framework  to  suppress  the  spurious  oscillation  in  the  vicinity  of 
sharp  thermal  fronts,  while  not  contaminating  the  smooth  flow  features  elsewhere.  The  resulting 
numerical  formulation  is  then  validated  on  several  benchmark  test  cases,  including  a  shock 
discontinuity  problem  with  the  ID  Burger’s  equation,  and  two  test  cases  for  the  compressible 
Euler  equations:  a  rising  thermal  bubble  and  density  current.  The  results  indicate  that  the 
proposed  DG-localized  Laplacian  artificial  viscosity  method  works  robustly  with  a  wide  range  of 
grid  sizes  and  polynomial  orders. 

1.  Introduction 

Numerical  weather  prediction  (NWP)  models  have  been  profoundly  influenced  by  the 
paradigm  shift  in  high  performance  computing  (HPC).  On  the  one  hand,  the  ever  increasing 
computing  power  allows  researchers  to  run  nonhydrostatic  (NH)  models  at  resolutions  finer  than 
10  km  [1];  on  the  other,  both  HPC  and  the  intrinsic  complex  physical  processes  in  NH  modeling 
pose  many  challenges  to  the  development  of  numerical  methods,  e.g.,  local  numerical  algorithms, 
high-order  accuracy,  geometric  flexibility,  etc.  The  discontinuous  Galerkin  (DG)  method  has 
been  proven  to  be  an  ideal  candidate  to  accommodate  these  challenges  [2].  One  example  is  the 
Nonhydrostatic  Unified  Model  of  the  Atmosphere  (NUMA)  [3,  4],  which  has  been  successfully 
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applied  to  3D  limited-area  modeling  on  distributed-memory  computers  with  a  large  number  of 
processors  [3,  4]  as  well  as  with  adaptive  mesh  refinement  (AMR)  in  2D  [5], 

Despite  the  success  in  NH  modeling  by  high-order  accurate  (i.e.,  order>2)  methods  [2,  6], 
robust  and  efficient  stabilization  of  sharp  flow  gradients  (e.g.,  thermal  fronts)  or  flow 
discontinuities  (e.g.,  shock)  remains  challenging  in  the  design  of  high-order  methods.  Arguably, 
the  two  most  frequently  adopted  methods  to  stabilize  the  high-order  methods  in  the  presence  of 
non-smooth  flow  features  are  limiters,  e.g.,  the  total  variation  bounded  (TVB)  limiter  in  the 
numerical  framework  of  Runge-Kutta  discontinuous  Galerkin  (RKDG)  [7],  and  artificial 
viscosity. 

In  the  numerical  simulation  of  nonhydrostatic  mesoscale  atmospheric  modeling,  very  high- 
order  polynomials  can  be  used  to  approximate  the  solution,  as  shown  in  Reference  [2].  Under 
this  scenario,  the  implementation  of  limiters  will  be  extremely  time-consuming.  Furthermore, 
after  limiting,  the  solution  might  be  represented  by  a  lower-order  or  even  piecewise  constant 
reconstruction.  This  polynomial  order  reduction  will  dramatically  increase  the  numerical 
dissipation  of  the  DG  algorithm  in  the  neighborhood  of  the  limited  element.  Sometimes,  key 
flow  features  can  be  totally  smeared  out,  especially  on  coarse  meshes.  Furthermore,  some  of  the 
most  effective  positivity-preserving  limiters  are  not  shape-preserving  [8].  Artificial  viscosity 
provides  an  alternative  way  to  handle  high-order  simulations  on  coarse  (i.e.,  under-resolved) 
meshes  in  the  presence  of  sharp  fronts. 

The  idea  of  capturing  shock  wave  discontinuities  in  a  fluid  by  adding  artificial  viscosity  into 
hyperbolic  conservation  laws  originated  from  Von  Neumann  and  Richtmyer  [9]  in  1950.  Since 
then,  many  types  of  artificial  viscosity  methods  have  been  developed  to  deal  with  flow 
discontinuity  capturing.  One  crucial  issue  in  all  artificial  viscosity  modeling  is  how  to  describe 
the  smoothness  of  the  flow  fields  accurately.  Smoothness  indicators  are  used  for  this  purpose. 
Different  smoothness  indicators  have  been  designed  based  on  the  gradient  of  flow  quantities  (e.g., 
velocity,  internal  energy,  etc.)  [10,  11],  the  resolution  of  numerical  representation  [12,  13],  the 
residual/entropy  residual  of  simulation  [14,  15,  16],  and  so  on.  Note  that  all  these  smoothness 
indicators  can  effectively  localize  the  artificial  viscosity  in  the  vicinity  of  flow  discontinuities. 
Based  on  the  different  procedures  to  design  artificial  diffusive  terms  and  to  incorporate  them  into 
the  original  governing  equations,  the  artificial  viscosity  methods  for  computational  fluid 
dynamics  can  be  roughly  classified  into  several  categories.  These  include,  but  are  not  limited  to 
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the  streamline-upwind/Petrov-Galerkin  (SUPG)  type  artificial  viscosity  [17,  18,  19,  20], 
localized  artificial  diffusivity  using  physical  principles  [10,  11,  21,  22,  23,  24,  25],  the  residual 
based  artificial  viscosity  [14,  15,  26,  27,  28],  the  entropy  artificial  viscosity  [16,  29,  30],  the 
spectral  vanishing  viscosity  [12,  31],  and  the  Laplacian  artificial  viscosity  [13,  32,  33,  34],  Other 
studies  of  the  artificial  viscosity  methods  can  be  found  in  References  [35,  36,  37,  38,  39],  just  to 
name  a  few. 

In  this  study,  considering  the  features  of  the  governing  equations  [2],  we  augment  the 
original  hyperbolic  system  with  the  localized  Laplacian  artificial  diffusive  terms  [13],  As 
mentioned  previously,  the  localized  Laplacian  artificial  viscosity  is  reconstructed  based  on  the 
smoothness  of  the  flow  fields.  Therefore,  an  adequate  amount  of  artificial  viscosity  is  localized 
in  the  vicinity  of  sharp  fronts  to  suppress  the  Gibbs  oscillation.  Meanwhile,  vanishing  artificial 
viscosity  does  not  contaminate  the  smooth  flow  features  away  from  sharp  fronts. 

The  paper  is  organized  as  follows.  The  governing  equations  for  the  nonhydrostatic  mesoscale 
atmospheric  modeling  and  the  discontinuous  Galerkin  discretization  are  introduced  in  Sec.  2.  In 
Sec.  3,  basic  ideas  behind  the  localized  Laplacian  artificial  viscosity  method  are  reviewed.  A 
new  family  of  modified  Laplacian  artificial  viscosity  models  is  introduced  based  on  the  proposed 
modeling  principles.  Sec.  4  then  presents  the  numerical  results  from  simulations  of  benchmark 
test  cases.  The  sensitivity  of  free  parameters  in  artificial  viscosity  modeling  is  also  studied  there. 
Finally,  conclusions  are  summarized  in  Sec.  5. 

2.  Governing  equations  and  discretization 

Many  different  forms  of  the  governing  equations  have  been  used  for  numerical  weather 
prediction  together  with  various  numerical  methods.  For  non-hydrostatic  atmospheric  modeling, 
three  equations  sets  were  presented  in  [2],  namely,  the  non-conservative  form  using  Exner 
pressure,  momentum,  and  potential  temperature  (Set  1),  the  conservative  form  using  density, 
momentum,  and  potential  temperature  (Set  2),  and  the  conservative  form  using  density, 
momentum,  and  total  energy  (Set  3).  It  was  found  in  [2]  that  the  two  conservative  forms 
outperform  the  non-conservative  form.  Therefore,  we  study  equation  Set  2  in  this  paper  which  is 
one  of  the  equation  sets  used  in  the  NUMA  model  [3,  4]. 

2.1  Governing  equations 
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(1) 


The  2D  form  of  equation  Set  2  reads 

dQ 

-^  +  V-F(Q)  =  Gm 

where  Q  =  (p,  pu,  pw,  p9)  are  the  conservative  variables,  p  is  the  density,  u  and  w  are  velocities 

in  X  and  z  directions,  respectively,  and  0  the  potential  temperature.  F  =  (/^,  /^)  is  the  inviscid 

flux  and  G  is  the  source  term.  They  are  defined  as 

pu  \  f  pw  \  10 

pv}  +  p 


yrx  _ 


puw 

pu9 


ir  = 


puw 
pw^  +  p 
pw6 


,  and  G  = 


(2) 


where  g  is  the  gravitational  constant,  p  is  the  pressure,  and  is  related  with  6  by  the  equation  of 
state  as  follows: 


p  =  Po 


(3) 


where  y  =  —  is  the  ratio  of  specific  heats  (for  constant  pressure 

Cp 


aim  (.;uii5LaiiL  vuiuiiic^. 


gas  constant,  and  po,  is  a  reference  pressure  that  is  only  a  function  of  the  vertical  coordinate. 
Introducing  the  splitting  of  the  density,  pressure  and  potential  temperature  as  p  =  po  +  p', 
p  =  Po  +  p',  and  9  =  9q  +  9',  where  the  subscript  ‘0’  denotes  the  values  in  hydrostatic  balance, 
we  rewrite  Eq.  (1)  as 

dO' 

-^+V-F'(Q)  =  G'(Q),  (4) 


where  Q'  =  (p',pu,pw,  &'),  6  =  p9  and  0'  =  0  —  po^o-  Correspondingly,  F'  is  written  as 


jrxf  _ 


,r  = 


and  G' 


(5) 


The  governing  equations  are  solved  on  the  physical  domain which  is  partitioned  into  N 
non-overlapping  elements The  solution  on  each  element belongs  toQ^(/2‘),  where 
is  the  space  of  tensor  product  of  polynomials  of  degree  at  most /c  in  each  variable 
defined  on/2'.  For  conciseness,  the  element-wise  continuous  solution  Q-  is  replaced  with  Q  in 
the  following  sections  when  no  confusion  between  Q-  and  Qi  is  present.  The  same  convention 


also  applies  to  F'  and  G'. 
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2.2  Discontinuous  Galerkin  method 

We  approximate  the  exact  solution  of  the  conservation  law  using  an  element-wise  continuous 
polynomial  G  =  [W  E  Let  W  be  an  arbitrary  weighting  function  or  test 

function  from  the  same  space  .  The  weighted  residual  form  of  the  governing  equations  on 
each  element  fl'-  then  reads 


Jo 


dt 


Jc 


WdV  +  V-F(Q^)fLdL 


Jn 


G(Qn)WdV, 


VW  E 


(6) 


(7) 


“'m  “'n‘ 

Applying  integration  by  parts  to  Eq.  (6),  one  obtains 

[  ^WdV-(  VW-F(Q„)dV+(  F-nWdS=(  G(Qi,)WdV, 

Jni  Jai  Jdiii  7ni 

where  F  =  (f.g)  and  n  is  the  outward  unit  normal  vector  of 

It  is  clear  that  the  surface  integral  in  Eq.  (7)  is  not  properly  defined  as  the  numerical  solution 

is  discontinuous  across  element  interfaces.  In  order  to  ensure  conservation,  the  normal  flux  term 

f  ■  n  is  replaced  with  a  Riemann  flux  F^jn(Ql,  Ql^,n),  where  denotes  the  solution  outside 

the  current  element  12'.  Various  (approximate)  Riemann  solvers  can  be  used  to  calculate  the 

Riemann  flux,  and  the  Rusanov  Riemann  solver  is  adopted  in  this  paper.  Then  Eq.  (7)  can  be 

rewritten  as 

“'n‘  “'5n‘  ■'n* 

In  the  DG  approach,  a  finite-dimensional  basis  set  [Wj]  is  chosen  as  the  solution  space.  Then 
the  governing  equation  is  projected  onto  each  member  of  the  basis  set.  Thus,  Eq.  (8)  is 
reformulated  as 
d 

>0}  ^  Jn,i  JdQ} 

(9) 


[  ^WdV-[  VW-F(Q^)dV+[  FJi^(QiQi;,n)WdS=  [  G(Q^)WdV. 

Joi  Joi  JdOi  Joi 


(8) 


y  (Qnjy^j)  dV  -  f  VlVfe  ■  F(Q^)dV  +  f  W,FJ^,^dS 
Jni  ^  JQ}  Jf)Oi 

=  [  W,y(GjWj)dV. 

Jni  ^ 


Applying  integration  by  parts  again  to  the  seeond  term  of  Eq.  (9),  the  strong  form  is  obtained 


as 


d 

dt 


f  ^ky(Qhji^j)dv+  [ 

Jn^  ^  Jn^ 


■  F(Qf^)dV  +  f  W,(FJi^ 


F^)dS 


(10) 
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=  I  ^  W„^{GjWj)dV, 

j 

where  F”  =  F  ■  n  is  the  local  flux  projected  on  dfl'  in  the  surface  normal  direction. 

The  first  integral  in  Eq.  (10)  is  usually  written  as  a  multiplication  of  the  mass  matrix  M  and 
the  time  derivative  of  the  solution  vector  [Q/j].  The  square  bracket  ‘[]’  denotes  the  vector  form  of 
the  solution  The  entries  of  the  mass  matrix  M  are  of  the  form 

WkWjdV.  (11) 

If  F  is  a  linear  function  of  Q,  then  F  can  be  expressed  as  F  =  Xy  FjWj.  Under  this  constraint,  the 
second  integral  in  Eq.  (10)  can  be  formulated  as  a  multiplication  of  the  stiffness  matrix  5*  and 
the  flux  vector  Ff  The  entries  of  the  stiffness  matrix  5*  are  written  as 

However,  if  F  is  a  nonlinear  function  of  Q,  then  F  cannot  generally  be  expressed  via  the  basis  set 
[Wj].  Quadratures  are  used  to  compute  the  volume  and  surface  integrals.  Clearly  these  operations 
can  be  expensive,  and  some  cost-effective  approaches  are  required  to  improve  the  computational 
efficiency.  One  such  solution  is  the  quadrature-free  approach  proposed  in  [40].  In  this  approach, 
it  is  assumed  that  the  flux  F  is  a  polynomial  which  belongs  to  the  same  space  Q^(/2Q  as  that  of 
the  solution  and  denote  it  by  F^.  Then  Eq.  (10)  still  holds  for  F^. 

We  also  assume  that  F^”^  belongs  to  the  polynomial  space  P^(df2Q  and  can  be  expressed  by 
the  basis  set  as  F^^y.  =  T.j  Fcom.fj^fj  surface.  Thus  mass  matrices  Bf  for  the 

surface  integration  in  Eq.  (10)  can  be  formed  with  entries 

=  f  (13) 

Substituting  Eqs.  (1 1)  -  (13)  into  Eq.  (10),  we  obtain  the  following  vector  form 

2 

^  -  Ff]  +  [Cj],  (14) 

1  =  1  f 

Now  consider  the  nodal  type  allocation  of  degrees  of  freedom  (DOFs)  [41],  and  assume  that 
Wjn  is  the  Lagrange  polynomial,  which  satisfies  14^ (ry)  =  S^j,  where  Vj  =  {xj,Zj)  is  the  nodal 
point.  Following  Ref.  [41],  we  introduce  the  differentiation  matrix  with  the  entries 
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(15) 


D 

x(0',m)  Q^l 

Then  the  entries  of  MD^i  can  be  calculated  as 


dWrr. 


dV 


dWr. 


^  dx^ 


dV 


(16) 


=  £  =  (s'W) 


Therefore,  Eq.  (14)  can  be  rewritten  as 


d[Qh\ 

dt 


=  “X  -  Ff]  +  [G>.]-  (17) 

1  =  1  f 


According  to  Eq.  (17),  in  the  implementation  of  the  strong  form,  there  is  no  need  to  explicitly 
calculate  the  stiffness  matrix  S’-,  but  the  differentiation  of  the  flux  polynomials.  This  fact  can  be 
utilized  to  save  computational  cost.  More  detailed  information  about  this  implementation  can  be 
found  in  Ref  [2]. 


3.  Localized  Laplacian  artificial  viscosity 

The  Laplacian  artificial  viscosity  is  used  to  suppress  the  Gibbs  oscillation  near  sharp  thermal 
fronts.  Generally,  for  2D  problems,  the  Laplacian  diffusion  terms  V  ■  F“’^((2,VQ)  in  x  and  z 
directions  read 


yrai?  _ 


du  \ 

du  \ 

dx 

^^■^P  dz 

dw 

and  = 

dw 

de'  i 
Ve,xP  } 

de’ 

Ve,zP  ) 

(18) 


For  simplicity,  we  set  Sg  x  =  Sg  ^  =  Sg. 

The  DG  method  is  used  to  discretize  the  following  equivalent  system  of  Eq.  (4), 

s-yQ  =  o, 

^  +  V  ■  F‘”^(Q)  -  V  ■  F^^{Q,S)  =  GdQ). 

Herein,  S  is  the  auxiliary  variable  to  facilitate  the  discretization  of  viscous  fluxes. 


(19) 
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The  artificial  viscosity  £  is  modeled  following  the  approach  in  Ref  [13],  Several 
modifications  are  introduced  to  make  this  model  more  suitable  for  sharp  thermal  front  capturing 
in  non-hydrostatic  atmospheric  modeling.  In  this  study,  the  resolution-based  indicator  is  used  to 
detect  non-smooth  flow  features.  Specifically,  we  approximate  the  solution  in  the  polynomial 


space  (f2)  as  follows. 


Q  «  ^  u^4>^. 


where  U  is  the  polynomial  approximation  of  Q,  (pi  is  the  ith  basis  of  the  space  and  N(k) 

is  the  total  number  of  basis  of  for  2D  problems,  N(k)  =  (k  +  l)x(k  +  1). 


Now  we  project  the  solution  U  onto  the  polynomial  space 

iV(fe-l) 


^(/2),  and  obtain 


=  ^  uSi. 


The  expansion  coefficients  Ui  can  be  calculated  by  solving  the  following  linear  system, 

iV(fc-l)  N{k) 

^  j  =  1,  -  ,  N(k  -  1)  . 

m=l  m=l 

Note  that  (■,■)  indicates  the  inner  product  in  L2(^2). 

The  resolution-based  indicator  in  one  finite  element  can  then  be  defined  as 


Sg 


{U  -UP,U  -UP), 
{U,U)g 


Finally,  a  smooth  variation  of  the  element-wise  artificial  viscosity  Eg  is  reconstructed  as 
follows, 

f  0  •£  c  ^  c  _ 

,  .  n(Sg-So)\  %  7  , 

fe  =  <y(l  +  sm - — - 1  if  So-K<Sg<So  +  K  (24) 

^  ^  ^  if  Sg>SQ  +  K. 

\  Eq 

It  is  clear  thatfg  G  [0,  £o].Note  that5o  in  Eq.  (24)  is  the  estimated  value  of  the  smoothness 
indicator  Sg  for  smooth  flow  features.  According  to  Ref.  [13],  if  the  polynomial  expansion  has  a 
similar  behavior  to  the  Fourier  expansion,  the  smoothness  indicator  will  be  proportional  to 
— 41o^io(^)-  Based  on  our  analyses,  this  estimation  can  add  unnecessary  numerical  dissipation 
to  relatively  smooth  flow  features.  Therefore,  Sq  is  set  as  —Slog^Q^k)  in  this  study.  The 
parameter  k  determines  the  smoothness  range  on  which  the  artificial  viscosity  functions. 
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Generally,  k  needs  to  be  chosen  sufficiently  large  so  as  to  ensure  a  sharp  front  capturing  with 
smooth  transition  to  flow  fields  nearby.  It  is  found  that  k  can  affect  the  performance  of  artificial 
viscosity  more  than  the  other  parameters  in  Eq.  (24)  do.  More  test  results  on  this  parameter  will 
be  discussed  in  the  following  section. 

Different  from  the  modeling  approach  presented  in  Ref  [13],  the  artificial  viscosity  Sq  is 
modeled  as  follows.  First  we  recall  the  definition  of  the  Peclet  number  Pe  for  a  diffusion 
process. 


LU 


Pe  =  — 


(25) 


a 

where  U  is  the  characteristic  speed,  L  the  characteristic  length,  and  a  the  diffusion  coefficient. 
The  artificial  viscosity  Sq  is  proportional  to  a.  In  Refs.  [13,  32],  U  is  set  as  the  maximum 
absolute  value  of  the  characteristic  speed  lAl^^x-  ^  is  the  sub-cell  grid  size  h/P,  where  h  is  the 
element  size,  and  P  is  the  polynomial  order,  fo  is  set  to  be  equivalent  to  a. 

In  this  work,  different  models  to  bridge  £o  and  a  are  proposed  to  make  the  modeling  of  the 
artificial  viscosity  Sq  less  sensitive  to  the  element  size  and  polynomial  order.  The  principles 
followed  in  this  approach  include: 

•  The  artificial  viscosity  £o  is  non-negative; 

•  When  the  resolution  of  the  numerical  scheme  is  infinite,  i.e.,  /i  ^  0  or  P  ^  oo,  the 
artificial  viscosity  £o  ^  0; 

•  The  modeling  is  compatible  with  the  classic  results  from  the  2"^^  order  accurate  (or 
equivalently  P^  reconstruction)  methods. 

Instead  of  using  the  uniform  assumption  of  the  sub-cell  grid  size  h/P,  we  redefine  the  length 
scale  in  Eq.  (25)  as  the  maximum  distance  between  two  adjacent  quadrature  points  in  the  element, 
which  is  written  as  /^h^ax  —  ^^max  '  K  where  A^max^  scaled  in  [0,1],  is  the  maximum  distance 
between  two  adjacent  quadrature  points  in  a  standard  ID  element.  Thus,  a  reads 

^^max  ,111 

^  ~  ^  '  l^lmax-  (26) 

A  general  model  for  the  artificial  viscosity  £o  can  then  be  written  as 

^0  ^  '  l^-lmax-  (22) 

We  now  focus  on  the  modeling  of  the  non-dimensional  function /(A^^a^.).  Following  Ref 
[37],  we  require  that  when  the  P^  reconstruction  is  used,  the  function  /  passes  the  point 
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(1,1/ Pe).  This  is  consistent  with  the  definition  of  a  for  the  2"^^  order  finite  volume  method.  Then 
we  show  one  way  to  determine  a  region  of  the  function  /  that  can  satisfy  the  proposed  modeling 
criteria.  It  is  observed  that  one  possible  upper  bound  of  the  function  /  can  be  written  as 


fi^^max)  -  Mmax  ^  [0-1]  (28) 

It  is  not  difficult  to  verify  that/(A^^ax)  >  0;  if  ^  0,  then  fo  ^  0;  andf(A^^ax) 

passes  the  point  (1,1/Pe).  One  possible  lower  bound  of  the  function  /  can  be  expressed  as 


/(A^max)  = 


fO,  0  <  A^rnax  <  1 
1 

^^max  —  1 


(29) 


This  region  is  shown  in  Figure  1  as  the  shadowed  area.  Note  that  the  linear  function 
/(A^max)  =  A^^aj;/Pe  recovers  the  choice  in  Ref  [13,  32].  Based  on  our  tests,  the  linear 


distribution 


=  (30) 

is  used  to  relate  Sq  with  a.  Finally,  the  artificial  viscosity  fo  is  defined  as 

We  note  that  the  artificial  viscosity  Eg  given  in  Eq.  (24)  is  an  element-wise  constant 
distribution.  It  is  obvious  that  Eg  has  a  jump  on  element  interfaces  if  the  element-wise  constant 
distribution  is  used.  For  quadrilateral  elements,  a  bilinear  distribution  can  be  constructed  by 
interpolating  the  four  vertex  artificial  viscosity  values  to  the  desired  quadrature  points.  The 
value  of  artificial  viscosity  on  a  specific  vertex  is  calculated  by  averaging  all  values  from  the 
neighboring  elements  which  share  the  vertex. 


4.  Results  and  discussions 

In  this  section,  we  test  the  localized  Laplacian  artificial  viscosity  method  using  several 
benchmark  problems  with  the  presence  of  shock  waves  or  sharp  thermal  fronts.  In  order  to 
evaluate  the  performance  of  artificial  viscosity  on  grids  with  different  resolution,  a  wide  range  of 
grid  sizes  and  polynomial  orders  is  tested  in  each  problem.  In  all  simulations,  Sq  in  Eq.  (24)  is 
selected  as  — 3Zo^io(^)  the  Peclet  number  Pe  is  fixed  at  2. 
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4.1  ID  Burger’s  equation  tests 

In  this  section,  we  test  the  efficacy  of  the  localized  artificial  viscosity  for  the  ID  Burger’s 
equation.  The  ID  inviscid  Burger’s  equation  augmented  by  an  artificial  diffusive  term  reads: 


dU  d  /I  \  d  (  dU\ 


(32) 


where  x  E  [-1.1].  Periodic  boundary  conditions  are  enforced  atx  =  —  1  andx  =  1.  The  initial 
conditions  are  defined  as  U(x,0)  =  Uq(x)  =  1  +  sin  (nx)/2.  According  to  Reference  [42],  a 
moving  shockwave  will  develop  after  t  =  2 /tt  under  the  given  initial  conditions.  An  element¬ 
wise  constant  distribution  of  Sg  is  used  to  stabilize  the  shock  wave.  In  all  simulations  presented 
in  this  section,  k  is  chosen  as  6. 

First  of  all,  the  results  of  different  artificial  viscosity  models  presented  in  Section  3  are 
compared.  The  results  are  shown  in  Figure  2.  In  Fig.  2,  ‘Log’  denotes  the  case  with  /(A^^ax)  = 
(1-  logA^max)/Pe;  ‘Linear(-)’  the  case  with =  (2  -  A^^ax)/^e;  ‘Constant’  the 
case  with  fil^^max)  —  l/^^l  and  ‘Linear(+)’  the  case  with  fil^^max)  —  ^^max/P^-  Simulations 
with  both  and  reconstructions  on  ten  elements  are  carried  out.  From  Fig.  2,  we  observe 
that  the  model  ‘Log’  is  the  most  dissipative  method  and  the  model  ‘Linear(+)’  is  the  least 
dissipative.  It  is  also  clear  that  the  performance  of  the  model  ‘Linear(+)’  is  sensitive  to  the 
polynomial  order,  while  that  of  the  other  models  is  not.  Based  on  this  observation,  the  model 
‘Linear(-)’  will  be  exclusively  used  in  all  simulations  in  the  rest  of  the  paper. 

Next  we  compare  the  results  with  P^  and  P^  reconstruction  on  different  grids.  The  solutions 
at  t  =  1  are  presented  in  Fig.  3.  The  corresponding  local  solution  errors  with  respect  to  the  exact 
solution  of  the  inviscid  Burger’s  equation  at  t  =  1  are  plotted  in  Fig.  4.  Several  observations  are 
summarized  as  follows.  From  Fig.  3,  we  find  that  the  localized  artificial  viscosity  works  robustly 
for  a  wide  range  of  high-order  reconstruction  (e.g.,  from  P^  to  P®  in  the  current  test).  For  all 
cases,  the  shock  is  captured  in  one  element.  From  Figs.  4  and  5,  it  is  clear  that  the  localized 
artificial  viscosity  does  not  contaminate  the  smooth  flow  features  away  from  the  shock,  but 
merely  concentrates  in  the  non-smooth  flow  regions  to  suppress  the  Gibbs  oscillation.  From  Fig. 
5,  we  observe  that  as  the  resolution  of  the  numerical  scheme  becomes  finer  (i.e.,  the  element  size 
becomes  smaller  or  the  order  of  the  reconstruction  polynomial  becomes  higher),  the  amount  of 
artificial  viscosity  localized  in  the  vicinity  of  the  shock  wave  becomes  smaller.  This  follows  the 
modeling  rules  as  stated  in  Sec.  3. 
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4.2  Rising  thermal  bubble 

The  rising  thermal  bubble  problem  is  driven  by  buoyancy  effects.  Specifically,  a  dry  warm 
bubble  rises  in  a  constant  potential  temperature  environment,  and  interacts  with  the  ambient  air 
during  this  process.  The  initial  potential  temperature  perturbation  is  given  as  follows  [2]: 


if  r>  Tc 

if  r  <  Tc 

(33) 

where  =  0.5°C,  =  250m,  r  =  and  (Xc,Zc)  =  (500, 300)  m  is  the 

initial  geometric  center  of  the  bubble.  The  hydrostatic  potential  temperature  6q  for  this  case  is 
300/C.  The  simulation  domain  is  (x,z)  G  [0,1000]^  m.  The  thermal  bubble  evolves  until 
t  =  700s.  Four  resolutions,  namely,  20m,  10m,  5m  and  3.5m,  as  presented  in  [2],  are  adopted  in 
the  simulations.  The  resolution  is  defined  as  L/ingyi^xk),  where  L  is  the  domain  size  in  the  x  or 
z  direction,  is  the  number  of  elements  in  the  corresponding  direction,  and  k  is  the 
polynomial  order.  Unless  explicitly  specified,  k  in  the  artificial  viscosity  model  is  set  as  0.5  in  all 
simulations  presented  in  this  section. 

4.2.1  Results  from  localized  artificial  viscosity 

The  maximum  and  minimum  potential  temperature  perturbations  O^ax  at  t  =  700s 

with  various  flow  field  resolutions  are  presented  in  Fig.  6.  Note  that  since  initially  6'  G  [0, 0.5], 
it  is  then  expected  that  during  the  evolution  of  the  thermal  bubble,  9'  is  bounded  in  this  range. 
From  Figure  6,  it  is  found  that  the  localized  Laplacian  artificial  viscosity  functions  perform  well 
for  a  wide  range  of  grid  sizes  and  polynomial  orders.  Only  small  overshoots  of  potential 
temperature  perturbation  show  up  in  the  results.  As  the  resolution  of  flow  fields  becomes  finer, 
the  numerical  dissipation  becomes  smaller.  Correspondingly,  both  maximum  and  minimum 
potential  temperature  perturbations  approach  the  theoretical  bounds. 

Then  the  effects  of  k  on  flow  field  features  are  studied  with  reconstruction  on  a  20x20 
mesh  (i.e.,  the  resolution  is  5m).  The  potential  temperature  perturbation  fields  with  different  k, 
namely,  0.5,  1,  2,  3,  and  4,  are  shown  in  Fig.  7.  It  is  observed  that  as  k  increases,  the  plume-like 
flow  features  near  the  thermal  front  are  gradually  damped.  From  the  maximum  and  minimum 
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potential  temperature  perturbations  0,^ax  at  t  =  700s  as  tabulated  in  Table  1,  it  is  clear 

that  the  overshoot  of  for  all  cases  is  very  small,  and  decreases  quickly  as  k  increases. 

The  mass  and  energy  conservation  properties  are  studied  for  low  resolution  cases,  including 
both  20m  and  10m  cases.  The  mass  and  energy  are  defined  as 


M(t)  =  j  p(t)dV  andE{t)=  f  p(t)e(t)dV, 


(34) 


where  e  is  the  total  energy.  In  this  case,  e  is  calculated  as  +  -  (u^  +  v^).  Correspondingly, 
the  mass  and  energy  loss  are  defined  as 


M_Loss(t:)  = 


M(t)  -  M(0) 


M(0) 


and  E_Loss(t)  = 


Eit)-Ei0) 


f(0) 


(35) 


The  results  for  solution  reconstruction  on  both  5x5  and  10x10  meshes  are  shown  in  Fig. 
8.  It  is  found  that  the  localized  artificial  viscosity  can  ensure  mass  conservation  and  only 
dissipates  internal  energy  which  is  to  be  expected  since  the  artificial  viscosity  used  here  is  not 
meant  to  represent  the  proper  Navier-Stokes  viscous  stress  terms. 


4.2.2  Comparison  between  localized  artificial  viscosity  and  limiters 

To  examine  the  advantage  of  the  localized  artificial  viscosity  method  on  handling  various 
high-order  simulations  on  coarse  meshes,  the  rising  thermal  bubble  case  is  run  with  low 
resolution  (i.e.,  20m  and  10m)  using  both  and  reconstruction.  The  results  are  then 
compared  with  those  from  a  limiter  using  the  combined  hierarchical  moment  limiting  procedure 
[43]  and  accuracy-preserving  positivity  limiting  procedure  [8]. 

A  minmod  TVB  (total-variation-bounded)  marker  based  on  the  potential  temperature  O'  is 
used  to  detect  the  “troubled”  cell  in  the  hierarchical  moment  limiting  procedure.  For  “troubled” 
quadrilateral  elements,  a  tensor  product  of  the  ID  mean-preserving  basis  [43]  is  used  to  carry  out 
the  solution  reconstruction.  The  maximum  polynomial  order  for  the  solution  reconstruction  in 
the  “troubled”  cells  is  fixed  at  two  (i.e.,  3’^'^  order  accurate).  In  the  accuracy-preserving  positivity 
limiting  process,  the  potential  temperature  O'  in  the  element  with  negative  O'  is  limited  as 
follows 

§'  =  a{0'  -  O')  +  O',  (36) 
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Qf  _ 

where  a  =  , — ,  6'  is  the  cell-averaged  value,  =  minj^sine}  is  set  of 

^  ~^min 

indices  of  all  quadrature  points  in  element  fig,  andf  is  a  small  positive  number  (e.g.,  10“^^  in 
this  study).  More  details  about  the  implementation  of  the  two  limiting  procedures  can  be  found  in 
[8,  43], 

The  potential  temperature  perturbation  fields  at  700s  from  simulations  using  localized 
artificial  viscosity  or  limiters  on  the  coarse  mesh  with  resolution  of  20m  are  displayed  in  Fig.  9. 

solution  reconstruction  is  used  for  (a)  and  (c)  with  localized  artificial  viscosity  and  limiters, 
respectively;  solution  reconstruction  is  used  for  (b)  and  (d)  with  localized  artificial  viscosity 
and  limiters,  respectively.  To  ensure  the  same  resolution  for  all  simulations,  a  17xl7mesh  is 
used  for  reconstruction,  and  5x5  for  reconstruction.  From  this  figure,  we  observe  that  the 
flow  fields  using  localized  artificial  viscosity  are  much  smoother  than  those  using  limiters.  The 
results  using  solution  reconstruction  with  P^  limiting  procedure  cannot  preserve  the  shape  of 
the  rising  thermal  bubble.  Similar  conclusions  can  be  drawn  from  Fig.  10,  which  shows  the 
potential  temperature  perturbation  fields  at  700s  with  similar  numerical  setup  as  that  in  Fig.  9, 
but  a  34x34  mesh  for  P^  reconstruction  and  10x10  for  P^°  reconstruction  (i.e.,  the  resolution  is 
10m).  All  these  results  demonstrate  the  superior  properties  of  localized  artificial  viscosity  on 
stabilizing  flows  with  thermal  fronts  for  a  wide  range  of  polynomial  orders  and  grid  sizes. 


4.2.3  Comparison  between  localized  artificial  viscosity  and  constant  viscosity 

Currently  a  common  practice  to  suppress  Gibbs  oscillation  in  thermal  front  capturing  is  to 
add  constant  viscosity  [44]  to  the  governing  equations.  Specifically,  the  physical  viscous 
diffusion  term  V  ■  F^(Q,S7Q)  is  added  to  the  right-hand  side  of  Eq.  (4).  F^(Q,S7Q)  in  x  and  z 
directions  can  be  written  as 


r  = 


/  0  \ 

du\ 

/  0  \ 

/  du\ 

PPTz 

dw 

and  g'^  = 

dw  , 

PP^ 

PPJI 

..30 

..30  / 

dx^ 


dz^ 


(37) 


where  p  is  the  constant  viscosity. 

It  is  obvious  that  this  approach  adds  numerical  dissipation  to  the  entire  flow  field,  no  matter 
whether  the  local  flow  features  are  smooth  or  not.  The  potential  temperature  perturbation  fields 


14 


at  700s  for  solution  reconstruction  on  a  10x10  mesh  using  a  series  of  constant  viscosity, 
namely,  O.lm^/s,  0.2m^/s,  0.3m^/s,  0.5m^/s,  Im^/s  and  2m^/s,  are  presented  in  Fig.  11.  The 
corresponding  maximum  and  minimum  potential  temperature  perturbations  ^'min 

t  =  700s  using  localized  artificial  viscosity  and  constant  viscosity  are  tabulated  in  Table  2.  From 
these  results,  we  observe  that  for  the  rising  thermal  bubble  case,  the  performance  of  constant 
viscosity  with  =  0.2m^/s  is  very  similar  to  that  of  localized  artificial  viscosity  as  shown  in  Fig. 
10(d).  If  the  constant  viscosity  is  very  large,  e.g.,  fi  =  2m^/s,  as  shown  in  Fig.  11(f),  the  flow 
structures  can  be  severely  dissipated.  Although  a  10x10  mesh  is  used,  the  resolution  of  the  case 
with  jU  =  2m^/s  is  very  similar  to  the  localized  artificial  viscosity  case  on  a  5x5  mesh  as  show 
in  Fig.  9(d).  More  advantages  of  the  localized  artificial  viscosity  approach  over  the  constant 
viscosity  approach  will  be  presented  in  Sec.  4.3.2. 

4.3  Density  current 

Now  we  study  the  density  current  problem.  In  this  case,  a  cold  bubble  drops  in  a  neutrally 
stratified  atmosphere,  hits  the  ground,  and  generates  Kelvin-Helmholtz  rotors.  The  initial 
potential  temperature  perturbation  is  given  as  follows  [2]: 

I  [l +  (38) 

where  =  — 15°C,  =  Im,  r  =  ^  (^c^^c)  =  (0,3000)  m  is  the  initial 

center  of  the  bubble,  and  =  (4000, 2000)  m.  Similarly  to  the  rising  thermal  bubble  case, 

the  hydrostatic  potential  temperature  Oq  is  set  to  300 K.  The  simulation  domain  is  (x,z)  G 
[0,2 5600]  X  [0,6400]  m.  The  cold  bubble  evolves  until  t  =  900s.  Four  resolutions,  namely, 
400m,  200m,  100m  and  50m,  are  used  in  the  simulations.  In  Ref  [2],  a  constant  dynamic 
viscosity  is  used  to  ensure  a  grid-converged  solution  at  approximately  50m  resolution.  Without 
explicit  viscosity,  the  simulation  will  eventually  blow  up.  We  now  present  a  flow  feature  based 
artificial  viscosity  to  stabilize  the  simulation.  Unless  explicitly  specified,  k  in  the  artificial 
viscosity  model  is  set  as  1  in  all  simulations  presented  in  this  section. 

4.3.1  Results  from  localized  artificial  viscosity 
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The  maximum  and  minimum  potential  temperature  perturbations  and  at  t  =  900s 
with  various  flow  field  resolutions  are  presented  in  Fig.  12.  Similar  conclusions  can  be  drawn 
from  this  figure  as  those  for  the  rising  thermal  bubble  case.  The  localized  Laplacian  artificial 
viscosity  works  well  in  a  wide  range  of  grid  sizes  and  polynomial  orders. 

The  effects  of  k  on  flow  field  features  are  studied  with  reconstruction  on  both  8x2  (i.e., 
400m  resolution)  and  64x16  meshes  (i.e.,  50m  resolution).  The  potential  temperature 
perturbation  fields  with  different  k,  namely,  0.25,  0.5,  and  1,  on  the  coarse  mesh,  and  those  with 
K  =  0.5, 1, 2, 4, 6  on  the  fine  mesh  are  displayed  in  Figs.  13  and  14,  respectively.  It  is  found  that 
the  artificial  viscosity  is  very  dissipative  on  the  coarse  mesh,  even  when  a  small  k  is  used.  For 
the  fine  grid  results,  as  k  increases,  fewer  Kelvin-Helmholtz  rotors  are  generated.  In  Table  2,  we 
tabulate  the  maximum  and  minimum  potential  temperature  perturbations  O^ax  ^'min 
t  =  900s  for  the  fine  grid  results.  It  is  clear  from  Table  2  that  the  overshoot  of  O^ax  for  cases 
is  small,  and  decreases  quickly  as  k  increases,  especially  when  k  exceeds  2. 

4.3.2  Comparison  between  localized  artificial  viscosity  and  constant  viscosity 

The  potential  temperature  perturbation  fields  at  900s  for  solution  reconstruction  on  a 
16x4  mesh  using  localized  artificial  viscosity  and  a  series  of  constant  viscosity,  namely, 
50m^/s,  75m^/s,  lOOm^/s  and  125m^/s,  are  displayed  in  Fig.  15.  The  corresponding 
maximum  and  minimum  potential  temperature  perturbations  O^^x  ^od  at  t  =  900s  using 
localized  artificial  viscosity  and  constant  viscosity  are  tabulated  in  Table  4.  A  similar  trend  can 
be  concluded  as  that  in  Sec.  4.2.3. 

It  is  found  that  if  the  constant  viscosity  is  “small”,  e.g.,  p  =  25m^/s,  the  simulation  diverged. 
Note  that  p  =  2m^/s  is  considered  a  “large”  viscosity  value  in  the  rising  thermal  bubble  case 
(This  is  true  even  if  we  consider  the  dimensionless  parameter  ^)-  Therefore,  from  the 

comparison  of  these  two  cases,  we  conclude  that  the  constant  viscosity  approach  suffers  from  the 
large  variation  of  viscosity  for  stabilization  purpose.  However,  the  value  of  localized  artificial 
viscosity  is  determined  by  the  numerical  resolution  of  the  scheme,  and  almost  no  parameter 
adjustment  is  needed  in  simulations  of  different  problems.  This  is  one  big  advantage  of  the 
localized  artificial  viscosity  approach  over  the  constant  viscosity  approach. 
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5.  Conclusions 

We  present  a  coupled  DG-localized  Laplacian  artificial  viscosity  method  to  suppress  Gibbs 
oscillation  near  sharp  thermal  fronts  in  nonhydrostatic  mesoscale  atmospheric  modeling. 
Specifically,  the  original  inviscid  governing  equations  are  augmented  by  Laplacian  artificial 
diffusive  terms.  The  diffusivity  is  a  function  of  the  local  smoothness  of  the  flow  fields.  Thus,  the 
proposed  method  has  a  favorable  sub-cell  shock  capturing  property,  and  does  not  contaminate 
the  smooth  flow  features  away  from  the  non-smooth  regions,  as  demonstrated  by  the  simulation 
results  for  the  ID  Burger’s  problem. 

In  order  to  alleviate  the  sensitivity  of  the  free  parameters  in  artificial  viscosity  modeling  on 
both  grid  sizes  and  polynomial  orders,  a  family  of  localized  artificial  viscosity  models  is 
proposed  and  tested.  We  use  this  numerical  framework  to  simulate  two  classical  2D  test  cases 
from  nonhydrostatic  mesoscale  atmospheric  modeling,  namely,  rising  thermal  bubble  and  density 
current  tests.  The  results  using  localized  artificial  viscosity  are  then  compared  with  those  using 
limiters  and  constant  viscosity.  The  results  show  that  the  proposed  artificial  viscosity  method 
works  robustly  with  a  wide  range  of  grid  sizes  and  polynomial  orders. 
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Figure  1.  Paradigm  of  the  family  of  functions  f(^^max)  Ihe  artificial  viscosity  model. 


(a)  (b) 

Figure  2.  Zoom-in  view  of  the  solutions  of  the  ID  Burger’s  equation  near  the  shock  wave  with  different 
artificial  viscosity  models  at  t  =  1  on  ten  elements,  (a)  reconstruction;  (b)  P^  reconstruction. 
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(a)  (b) 

Figure  3.  Solutions  of  the  ID  Burger’s  equation  at  t  =  1  on  different  grids,  (a)  reconstruction;  (b) 
reconstruction. 


(a)  (b) 

Figure  4.  Local  error  of  computed  solutions  of  the  ID  Burger’s  equation  at  t  =  1  on  different  grids,  (a) 
P^  reconstruction;  (b)  P^  reconstruction. 
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Figure  5.  Distribution  of  the  artificial  viscosity  from  the  ID  Burger’s  equation  simulation  at  t  =  1  on 
different  grids,  (a)  reconstruction;  (b)  reconstruction. 


(a)  (b) 

Figure  6.  The  maximum  and  minimum  potential  temperature  perturbations  0max  ^min  ^he  rising 
thermal  bubble  att  =  700s  with  various  flow  field  resolutions,  (a)  Blnax  degree  of  polynomial;  (b) 
^min  degree  of  polynomial. 
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Figure  7.  Potential  temperature  perturbation  fields  of  the  rising  thermal  bubble  at  t  =  7005  for  different 
K  with  reconstruction  on  the  20x20  mesh. 
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Table  1.  The  maximum  and  minimum  potential  temperature  perturbations  Oly^ax  ^min  ^^e  rising 
thermal  bubble  ait  =  700s  for  different  k  with  reconstruction  on  a  20x20  mesh. 


(a)  (b) 

Figure  8.  Conservation  of  (a)  mass  and  (b)  energy  for  the  rising  thermal  bubble  simulations  using 
localized  artificial  viscosity  on  two  different  meshes  with  P^^  reconstruction. 
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Figure  9.  Potential  temperature  perturbation  fields  of  the  rising  thermal  bubble  at  t  =  7005  using 
localized  artificial  viscosity  and  limiters  with  20m  resolution,  (a)  reconstruction  on  a  17x17  mesh 
with  limiters;  (b)  reconstruction  on  a  5x5  mesh  with  limiters;  (c)  P^  reconstruction  on  a  17x17 
mesh  with  localized  artificial  viscosity;  (b)  P^^  reconstruction  on  a  5x5  mesh  with  localized  artificial 
viscosity. 
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localized  artificial  viscosity  and  limiters  with  10m  resolution,  (a)  reconstruction  on  a  34x34  mesh 
with  limiters;  (b)  reconstruction  on  a  10x10  mesh  with  limiters;  (c)  P^  reconstruction  on  a  34x34 
mesh  with  localized  artificial  viscosity;  (b)  P^^  reconstruction  on  a  10x10  mesh  with  localized  artificial 


viscosity. 
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Figure  11.  Potential  temperature  perturbation  fields  of  the  rising  thermal  bubble  at  t  =  700s  for  constant 
viscosity  with  different  ji  using  reconstruction  on  the  10x10  mesh. 
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Table  2.  The  maximum  and  minimum  potential  temperature  perturbations  Oly^ax  ^min  ^^e  rising 
thermal  bubble  at  t  =  700s  for  localized  artificial  viscosity  and  constant  viscosity  with  different  ji  using 
reconstruction  on  a  10x10  mesh.  LAV  stands  for  localized  artificial  viscosity. 
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(a)  (b) 

Figure  12.  The  maximum  and  minimum  potential  temperature  perturbations  and  of  the 
density  eurrent  flow  at  t  =  900s  with  various  flow  field  resolutions,  (a)  9^ax  vs.  degree  of  polynomial; 
(b)  dlnin  vs.  degree  of  polynomial. 
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(c)  K  =  1 

Figure  13.  Potential  temperature  perturbation  fields  of  the  density  current  at  t  =  9005  for  different  k 
with  reconstruction  on  the  8x2  mesh. 
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(q)  K  =  6 

Figure  14.  Potential  temperature  perturbation  fields  of  the  density  current  at  t  =  900s  for  different  k 
with  reconstruction  on  the  64x16  mesh. 
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Table  3.  The  maximum  and  minimum  potential  temperature  perturbations  9max  ^min  ^^e  density 
current  at  t  =  900s  for  different  k  with  reconstruction  on  a  64x16  mesh. 
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Figure  15.  Potential  temperature  perturbation  fields  of  the  rising  thermal  bubble  flow  at  t  =  7005  for 
localized  artificial  viscosity  and  constant  viscosity  with  different  {i  using  reconstruction  on  the  16x4 
mesh.  LAV  in  (a)  stands  for  localized  artificial  viscosity. 
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9.243x10-2 

-8.835 

Table  4.  The  maximum  and  minimum  potential  temperature  perturbations  9max  ^min  ^^e  density 
current  at  t  =  9005  for  localized  artificial  viscosity  and  constant  viscosity  with  different  ji  using 
reconstruction  on  a  16x4  mesh.  LAV  stands  for  localized  artificial  viscosity. 
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